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ABSTRACT 

We explore the dependence of the subhalo mass function on the spectral index n of the lin- 
ear matter power spectrum using scale-free Einstein-de Sitter simulations with n = — 1 and 
n = —2.5. We carefully consider finite volume effects that may call into question previous 
simulations of n < —2 power spectra. Subhaloes are found using a 6D friends-of-friends al- 
gorithm in all haloes originating from high-cr peaks. For n = — 1, we find that the cumulative 
subhalo mass function is independent of the parameters used in the subhalo finding algorithm 
and is consistent with the subhalo mass function found in ACDM simulations. In particular, 
the subhalo mass function is well fit by a power-law with an index of a = —0.9, that is the 
mass function has roughly equal mass in subhaloes per logarithmic interval in subhalo mass. 
Conversely, for n — —2.5, the algorithm parameters affect the subhalo mass function since 
subhaloes are more triaxial with less well defined boundaries. We find that the index a is 
generally larger with a > —0.75. We infer that although the subhalo mass function appears 
to be independent of n so long as n > —2, it begins to flatten as n — > —3. Thus, the common 
practice of using a as —1.0 may greatly overestimate the number of subhaloes at the smallest 
scales in the CDM hierarchy. 

Key words: methods: iV-body simulations - methods: numerical - galaxies: haloes - galax- 
ies: subhaloes - dark matter 



1 INTRODUCTION 

In the current Cold Dark Matter (CDM) paradigm, structure forms 
hierarchically; small-scale density fluctuations collapse to form an 
early generation of haloes, which are the progenitors of larger-scale 
systems. While some progenitors survive as substructure, others are 
tidally disrupted and become the smooth component of the new sys- 
tem. The distribution of substructure within galaxy-size haloes is of 
great interest as galaxies and satellites, acting as baryonic tracers of 
the underlying mass distribution, can provide observational checks 
of the current CDM paradigm. The subhalo distribution at much 
smaller scales, likely devoid of baryonic tracers, is of practical in- 
terest for DM detection experiments. 

In this paper, we examine the properties of subhaloes in a pair 
of scale-free cosmologies meant to mimic two different scales in 
the Universe. Objects in a scale-free simulation can be related to 
objects at different scales in the CDM hierarchy via the scale de- 
pendence of effective spectral index n e fl = d In P(k) /d In k of the 
ACDM concordance model. We use scale-free simulations for their 
simplicity as there is only one physical scale in the simulation, the 
nonlinear scale. This feature allows us to examine whether halo 
substructure remembers initial conditions, namely the index of the 
initial power spectrum. 

Dark matter haloes in ACDM simulations exhibit properties 
of self-similar and fractal systems. For example galactic haloes ap- 



pear to be rescaled versions of cluster haloes JMoore et al.|1999] >. 
Moreover haloes contain subhaloes whose properties such as the 
density profile are similar to those of haloes (e.g. |Diemand et al.| 
2008, and Sprin gel et al.|2008) . The subhalo mass and circular ve- 
locity distributions show little scale dependence if normalized in 
terms of the host halo mass or maximum velocity (e.g. |Kravtsov] 
|et al.||2004||Gao et al.|2004| and |Reed et al.|2005| >. Furthermore, 
the subhalo mass distribution is characterized by a power-law, 
dN/dlnM sub oc Mf ub with a = -0.8 to -1.0 (e.g. |Stoehr 



[eTaLl[2003l |Gao et al.||2004l |Diemand et aTJ|2007l |Madau et al 
|2008| and |Springel et a l. 2008 1. Note that a power-law index of -1 
implies a scale-free distribution with a constant mass in substruc- 
ture per logarithmic mass interval. The most recent studies seem to 
favour a — —0.9 for galactic hal oes down to a subha lo mass of 



10 5 M jMadau et al. 



2008 



and 



Springel et al.|2008) . Substruc- 



ture is not entirely independent of environment as the amount of it 
in a given halo depends on the halo's formation time or peak height 
and concentration; haloes that form earlier from high-cr peaks have 
more substructure than 1ow-<t haloes of similar mass while the more 
concentrated the halo's density profile the more substructure it con- 
tains (e.g. |Bullock et al.|200f]|Gao et al.|2004| and |Zentner et al.| 
[2005> . 

These ideas suggest a simple picture of substructure: haloes 
contain a scale-free distribution of subhaloes, which in turn contain 
a similar distribution of subsubhaloes, all the way down the CDM 
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hierarchy. The scale of the bottom of the hierarchy is set by funda- 
mental physics, namely the free-streaming and collisional dampen- 
ing scales of dark matter. If, for example, dark matter is the neu- 
tralino, a Weakly Interacting Massive Particle (WIMP) predicted 
by supersymmetric extensions to the Standard Model (SUSY), 
the first objects form at a redshift of z > 60 and have masses 



M < 10" 6 M ( jGreen et al.|2004[|2005p 



Diemand et al. 



{2005} 



investigated the formation of these objects and postulated that there 
would be ~ 10 15 of them in the Milky Way (MW) halo today. The 
formation of a rare, high redshift 0.014 Mq object was simulated 
by Diemand et al. (2006). They found that the subhalo mass func- 
tion is very similar to that of cluster haloes. Again the suggestion is 
that self-similarity in halo's substructure extends all the way down 
the hierarchy. 

Self-similarity in the distribution of substructure may have im- 
portant implications for a variety of experiments such as GLAST, 
which will search for 7-rays from DM self-annihilation. The an- 
nihilation signal is very sensitive to the slope of the subhalo mass 
distribution and the number of WIMP-scale subhaloes since sub- 



structure can boost the flux by factors of a 3-100 (e.g. Diemand 



et al. 



2007)|Kuhlen et al.|2008[|Pieri et al.|2008"l and |Strigari etal 



2008 ). The amount of dark matter locked up in substructure also 
has ramifications for direct DM detection as it reduces the den- 
sity of the smooth background component of the Galactic halo but 
also increases the local density inside subhaloes ( Kamionko wski &] 
Koushiappas 2008 ). 

The extrapolation used to make predictions of the 7-ray flux 
are non-trivial since the MW halo is a factor of 10 18 times more 
massive than the smallest subhaloes. The estimate by |Diemand] 



et al 



(2005) of ~ 10 small subhaloes neglects mergers and 



tidal interactions inherent in the hierarchical structure formation 
scenario. Haloes must survive similar-mass mergers and accre- 
tion, along with other dynamical processes that exist in galaxies 
l |Zhao et aTj 2007 ). Furthermore, the choice of a is crucial as small 
changes in a by as little as 0.1 can change the number of subhaloes 
at the bottom of the hierarchy by an order of magnitude or more. 

These results suggest that the internal properties of individual 
haloes (e.g. the subhalo mass function) have no memory of the ini- 
tial power spectrum. However, this would be suprising considering 
other internal properties of dark matter haloes, such as concentra- 
tion of the density profile, depend on properties of the primordial 
power spectrum (e.g. Reed et aT|2005| >. For example, scale-free 
simulations by Knollmann et al. ( 2008 ) show that haloes have less 
concentrated density profiles as n — > —3. 

There are other theoretical reasons to suspect that substruc- 
ture is qualitatively different at small scales. As one approaches 
the bottom of the CDM hierarchy, n e ff monotonically decreases to 
—3 down to the cutoff at the WIMP free-streaming scale. The di- 
mensionless power spectrum, A 2 (fc) oc fc ncrr+3 , becomes scale in- 
dependent as n — > —3 and objects collapse simultaneously over a 
wide range in scales above the free-streaming scale. Haloes at these 
scales do not form in a clean hierarchical fashion and may not viri- 
alize before merging with or being accreted by other haloes, which 
may influence the distribution of substructure. In short, structure 
formation at the smallest scales in the CDM hierarchy may be qual- 
itatively different from structure formation via hierarchical cluster- 
ing that occurs at galactic scales. 

There are hints in previous studies that substructure does have 
a spectral dependence. Die mand et al.| ( (2006| l showed substructure 
in a high redshift subsolar mass object was more susceptible to tidal 
disruption, with as much as 20—40% of it being disrupted within an 
expansion factor of 1.3 as compared to 1% in a low redshift cluster 



halo. |Springel et al.| ([2008 1 found that subgalactic subhaloes have 
less substructure than galactic haloes when regions of the same 
overdensity are compared. They also presented evidence that sug- 
gests the slope of the subsubhalo mass function flattens as the sub- 
halo host mass decreases. The scale-free simulations of Re ed et aTl 
(2005 ) appear to show that the subhalo velocity distribution flattens 
when n < —2. 

However, simulations of n < —2 cosmologies are notoriously 
difficult to perform due to finite volume effects ( |Smith et al.|2003| >. 
Missing power from modes larger than the simulation box affect 
statistical quantities such as the two-point correlation function and 
halo mass function (Bagla & Ray 2005 , Power & Knebe 2006, and 
Bagla & Prasad 2006). Provided large-scale modes have negligi- 
ble amplitudes and the scale of interest is a small fraction of the 
box size, the internal properties of haloes appear unaffected. One 
can also correct for deviations in statistical quantities, though these 
corrections become increasingly large as n — > —3. Generally, it 
is preferable to simply halt a simulation before finite volume ef- 
fects become an issue. We revisit the criterion presented by [Smith] 
|etaT]{2003"l > that ensures the finite volume effects are negligible and 
show that it has not been met in the prior simulations of scale-free 
power spectra with n < — 2. 

We run a pair of scale-free Einstein-de Sitter simulations in 
order to quantify the spectral dependence of the subhalo mass 
function. As most studies focus on clusters and galaxies, where 
n e ff « — 1.8 and —2.1 respectively, we choose spectral indices of 
n = — 1 and n = —2.5 to bracket these scales in the CDM hierar- 
chy. Finite volume effects are carefully considered in choosing the 
end point of our simulations. We search all large haloes for sub- 
haloes and vary the parameters of the group finding algorithm in 
order to quantify systematics. Our results show that subhalo mass 
function does depend on n, in contrast to previous results, and is 
sensitive to the parameters used to find subhaloes. 

Our paper is organized as follows: In section we discuss 
particle number requirements and the difficulties with simulating 
n — > —3 cosmologies. The initial conditions and cosmological pa- 
rameters used along with the technical details of the simulations 
are presented in [[3] In ^4] we compare the halo mass distribution 
from our simulations with theoretical models. In Sj5] we present 
our results for the properties of subhaloes and the dependence of 
subhaloes on the spectral index. The paper concludes in Sj6]with a 
summary and discussion. 



2 PARTICLE REQUIREMENTS 

Simulations are a compromise between two competing goals: good 
statistics, which favours the use of a physically large box, and high 
resolution, which promotes the use of a small box to achieve high 
resolution of individual objects. Since substructure is sensitive to 
the softening length used in simulations the second of these issues 
tends to dominate the choice of box size. Indeed, substructure in 
haloes composed of < 10 particles tends to evaporates due to nu- 
merical softening effects (e.g. Moore et al.|1999| and |Klypin et al.| 
[T999l >. 

Choosing between a small box with a resolution high enough 
to resolve small-scale structure while still keeping large-scale 
modes in the linear regime is problematic. Simulations must be 
halted before the large-scale modes are nonlinear so as to ensure 
that mode-mode coupling is correctly represented, and that the 
amount of power due to missing large-scale modes is negligible. 
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Figure 1. Number of particles required, versus spectral index n 

based on Eq. (top) and Eq. \10) (bottom) for a variety of constraints 
on the linearity large-scale modes and the nonlinearity of modes at the 
Nyquist scale. The constraints for the top panel are (A 2 (fcNy), f^ iss ) = 
[(2,0.04), (1.0,0.1), (0.5,0.2)] in solid, dashed and dotted curves respec- 
tively. The bottom panel indicates the requirements for haloes of 10 5 par- 
ticles to be lcr, 2ct and 3<r peaks in solid, dashed and dotted curves re- 
spectively with the box mass being a 5cr peak. Filled stars indicate the 
resolution and spectral index of the two simulations we ran. For compar- 
ison, also shown are scale-free simulations from Knollmann et al. 1 2008 1 in 
open crosses, Reed et al. (2005 ) in open squares, and Jain & Bertschinger 
(1998) in open triangles. The effective resolution and effective spectral in- 
dex of the CDM power spect rum from the Nyqu ist scale up to fchaio = 



Diemand et al. 



(2006 



(47rp bg /3M) 1 /3 for ^ 
thick horizontal lines respectively. 



The dimensionless power spectrum is 



20071 



are shown in thin and 



A 2 (fc) = 



V 



(1) 



where V is the normalization volume and P(k) is the power spec- 
trum. In what follows, P(k) and A 2 (fc) refer to the power spectrum 
evolved according to linear theory. We assume a scale-free initial 



P(k) = a 2 Ak n , 



(2) 



where A is the amplitude and a is the cosmological scale factor. 
Modes are considered to be linear if A 2 (fc) <C 1 and nonlinear 
if A 2 (fc) ^ 1 with the nonlinear scale defined by the relation 
A 2 (A;nl) = 1- The effective index of the power spectrum is 



n eff (fc) 



rim A 
dink 



(3) 



For a simulation of size L with N A particles, the power at a given k 
is related to the power at the Nyquist wavenumber, &Ny = irN/ L, 
by 



A 2 (fc) = A 2 (fc Ny ) 



Ny 



3+n 



Thus, the power at the box scale kh = 2n/L is 



A 2 (fc b ) 



A 2 (fc Ny ; 



N 



3+n 



(4) 



(5) 



We can use this definition to define an end point for a given simu- 
lation by requiring that the mode at the box scale is still linear. 

The impact of modes larger than the box can be quantified 
by considering missing variance, that is the difference between the 
variance integral for an infinite universe and the finite sum over 
modes within the box. |Smith et a l. (2003 ) show that the missing 
variance for scale-free power spectra is well approximated by 



A 2 (k b ) 
3 + n 



F{3 + n) 



(6) 



where F(x) = 1 - 0.31s + 0.015x 2 + 0.00133a; 3 for -3 sC 
n ^ 1. As n — > — 3, <r 2 iss — > oo, that is, missing power plays an 
increasingly important role. Ideally, one wants cr, 2 iss — > 0, though 
so long as er 2 ^ <C 1 the simulation will not suffer significantly 
from finite volume effects and the large-scale modes will still be 
linear. 



Smith et al. 



(2003 



use (T^iss < 0.04 in their study of the 
nonlinear evolution of the power spectrum, which we also use as 
our strictest limit. However, there is not a precise value of (j 2 iss that 
guarantees some level of accuracy will be achieved for any given 
realization. 

We can combine a desired nonlinearity at a given scale and 
our constraint on <r 2 iss to obtain a minimum requirement on the 
number of particles used in a simulation. The minimum number of 
particles required to achieve a level of nonlinearity at the Nyquist 
scale, A 2 (fcN y ) is 



A 2 (fc Ny ) F(3 + n) 
^ 3 + n 



3/(3+n) 



(7) 



A useful choice is A (&Ny) > 1 since this ensures that the sim- 
ulation is in the nonlinear regime. It is evident from this equation 
that for A 2 (A;N y ) / <T 2 iss > 1, the minimum number of particles rises 
super-exponentially asrn —3. 

Another measure of nonlinearity is the mass variance, a 2 (M), 
given by 



\M) = 



V 



(27T) 



P{k)\W 2 (kR)\d 3 k, 



(8) 



where W{x) = (3/x 3 )(sin x — x cos x) is the Fourier transform 
of the top-hat window function, R = (3M/47rpb g ) 1 ^ 3 and pbg is 
the background mass density. A mass scale is linear if <r(M) < <5 SC , 
where <5 SC is the critical density for spherical collapse, and nonlinear 
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if a(M) > Ssc- The characteristic mass scale M* is defined by the 
relation c(M») = 5 K . We require <r(Mb x) <C &c, where Mb 0X = 
?rA^ 3 /6 is the mass enclosed in the largest sphere contain in the 
simulation volume. As a(M) is a statistical quantity, the criteria 
that (j(Mbox) = S^/u, where where v is the height of density field 
in rms units, is a probabilistic one. For example, with v — 2 and 
Gaussian initial conditions, 4.5% of an ensemble of simulations 
would be nonlinear at the box scale. We choose Vb m ~ 5, ensuring 
that the probability that the box mode is nonlinear is ~ 10~ 6 . Note 
that the effective index at a given mass scale is 



n cff (M) 



rilna 2 (M) 
d\nM 



(9) 



Using the mass variance, one can impose evolutionary criteria 
to calculate the minimum number of particles required. Generally, 
studies focus on haloes of a particular mass scale, M, which cor- 
respond to density peaks above some threshold v. The minimum 
number of particles required to investigate haloes originating from 
VhO peaks of composed of Nh particles is: 



AT 3 =^ 

1 y req 

IT 



'-'box 
Vh. 



6/(3+n) 



(10) 



M[wy 




k [Mpc ] 



As n — > —3, the mass variance becomes independent of scale and 
the minimum number of particles again rises super-exponentially. 

We demonstrate the relative importance of these constraints 
in Fig. [T] by plotting A?" r 3 q as a function of n. The curves indicate 
the minimum number required for a given set of constraints, with 
the constraints becoming relaxed as one goes from solid to dashed 
to dotted. To be representative of a given cosmology and spec- 
tral index, simulations must lie above these curves. Note that the 
number of particles needed to attain a highly evolved simulation 
A 2 (MMy) > 1 while still limiting o- 2 iss < 0.04 for n < -2 be- 
comes increasingly impractical. For very negative indices, a more 
reasonable goal is evolving the simulation till the Nyquist scale is 
just nonlinear while <7 2 is! . ~ 0.10. Again, we stress that though our 
choices of (T^ iss and Ubo X are meant to be conservative, these crite- 
ria do not guarantee that our study is free of finite volume effects. 
Simulations are always missing power and this missing power, even 
if small, can subtly affect the evolution of the system as shown in 
Widrow et al. (2009). A rigorous examination of finite volume ef- 
fects on the subhalo population will require further study. 

The missing variance is not as large in ACDM as it is for 
scale-free power spectra due to the turnover in the power spec- 
trum but still goes as A 2 (fcb). This figure clearly shows the dif- 
ficulties with modelling the bottom of the CDM hierarchy where 
ricff fa —2.8. The highest resolution simulation of SUSY haloes by 
|Diemand~et al, (2006)has n e ff ~ —2.85 to —2.75. As they examine 
the substructure of a very rare 3.5a peak and halt their simulation 
early enough, the missing variance does not exceed « 0.05. How- 
ever, any extrapolations based on a 3.5cr peak should be treated 
with caution. 

Equations (J7]l and < | 1 0[ > can also be rearranged to determine 
when to halt a simulation for a given <r^ iss or <r 2 (Mbox). Since we 
are interested in substructure, we focus on > 10 5 particles for n = 
-2.5 and N 3 = 720 3 , halting the simulation when o(M hm ) « 
<5 sc /4.9 means that these haloes originate from rare ~ 2.5<r peaks. 
If our n = —1 simulation is evolved to the same cr(A'/b ox ), the 
haloes composed of 10 5 particles would originate from common ~ 
0.4cr peaks. Using a 2 (M) oc a 2 , it is a simple matter to calculate 
how many e-foldings the simulation must be evolved to achieve this 
desired end state. 



Figure 2. Effective spectral index of the CDM power spectrum versus 
mass (solid) and wavenumber (dashed). Highlighted in thick lines are the 
wavenumbers and mass scales sampled in simulations such as the Via 
Lactea simulation and the Aquarius Project. Various mass scales in the 
CDM hierarchy are indicated by filled circles, with the largest mass corre- 
sponds to the an MW halo (10 12 Mq) and the other two points indicating 
the largest and smallest galactic subhalo masses examined in previous stud- 
ies, 10 10 Mq and 10 6 Mq respectively. Also shown are the indices of our 
n = — 1 and n = —2.5 simulations (dotted) with the filled circles indicat- 
ing masses of 5 X 10 5 , 10 4 and 100 particles corresponding to the mass of 
haloes , largest subhaloes and smallest subhaloes examined. 



3 NUMERICAL METHODS 

3.1 How to interpret scale-free simulations 

The varying slope of the ACDM power spectrum is key to inter- 
preting scale-free simulations in the context of the CDM hierar- 
chy. In Fig. [2] we show the effective spectral index as a function 
of wavenumber and mass based on the ACDM power spectrum 
of |Hu & Eisenstein ( 1998). Highlighted are scales sampled in the 
"Via Lactea" (VL) simulation along with certain mass scales. Most 
numerical studies focus on clusters and galaxies, corresponding to 
Ties ~ —1.8 to —2.2. Also shown are the indices of our simula- 
tions, where the scale of the box size is set to match the scale at 
which the CDM power spectrum has the same index. Our choice of 
indices is meant to bracket galaxy and cluster scales. In essence, 
haloes found in an n — —2.5 simulation are representative of 
high redshift haloes surrounding the first protogalaxies. Objects at 
this scale may also become the larger subhaloes found in galaxies. 
Though the actual power spectrum is steeper at these scales with 
n e ff w —2.7, we satisfy ourselves with n — —2.5 as it is the steep- 
est spectra we can reasonably simulate, though even at this index 
we are pushing the limits outlined in ^2] 



3.2 Simulations 

We run scale-free Einstein-de Sitter (fi m = 1) cosmological sim- 
ulations with power-law power spectra as shown in Eq. |2j. We 
choose an amplitude and initial scale factor so that the maximum 
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Table 1. Summary of Simulations Run. 



n TV 3 (af/m) a i Softening Length e s 



-1 

-2.5 



720 3 
720 3 



42.8 
17.5 



0.214 
0.564 



1/30(L/N) 
1/30(L/N) 



initial displacement for any given particle is less than 1/2 the ini- 
tial inter-particle spacing. Due to the random nature of the initial 
density field this choice does not correspond precisely to a specific 
amplitude requirement across all power spectra, but ensures that 
the simulation starts in the linear regime. We draw both simula- 
tions from the same random realization and normalize the ampli- 
tudes so that a(M = 15 x 10 6 particles, a = 1.0) = 0.9, which 
is equivalent to normalizing using as- A final issue concerns tran- 
sients in the density and velocity field generated by initial condi- 
tions generator. To minimize the effect of transients, we use ini- 
tial conditions generated by second-order Lagrangian perturbation 
theory (2LPT) rather than the standard Zeldovich approximation 
(ZA) (JCrocce et al. 2006). The simulations are run using the paral- 
lel N-body tree-PM code GADGET-2 l [Springel|2005l >. We halt the 
n — —2.5 simulation when haloes composed of > 10° particles 
originate from > 2.5<r peaks and the box scale corresponds to a rare 
4.9cr. Evolving farther in order to form larger haloes and further re- 
duce numerical softening effects is not possible since, even at the 
our chosen end point, we are pushing the limits with a^iss ~ 0.10. 
We also analyze the n = — 1 simulation when > 2.5<r haloes are 
composed of > 10 5 particles for the sake of consistency, though 



°miss ~ 0.001 at this point. A summary of the simulations is given 
in Table[[] 



4 HALOES 

We identify haloes at each time step using a friends-of-friends 
(FOF) group finder with a linking length ih — 0.2 times the co- 



moving interparticle spacing Ax = L/N (Davis et al. 1985 1. Only 
FOF groups with more than 32 particles are kept in the halo cat- 
alogue. The comoving halo number density is shown in Fig. [3] 
at three different redshifts along with theoretical predictions from 
|Sheth & Tomien1p999t (ST) and |Press & Schechter|fi974l > (PS). 



The 



PS 



mass function is based on the spherical collapse model 
while the |ST| incorporates ellipsoidal collapse with additional mass 
function parameters calibrated using large-scale ACDM simula- 
tions. For both simulations, the |ST| mass distribution is in better 
agreement than|PS|for all redshifts. 



5 SUBSTRUCTURE 

5.1 Searching for Subhaloes 

Identifying subhaloes is more difficult than identifying haloes. The 
key issue with defining the outer boundary of a subhalo embedded 
in a gravitationally bound object. Several methods have been devel- 
oped to search for substructure, including SKID (Stadel 2001 1 and 
SUBFIND (Springel et al. 2001). In this paper, we use a phase- 
space FOF (6DFOF) algorithm, which is similar to the one de- 
scribed in Diem and et al.| l [2006| l. We limit our subhalo search to 
haloes originating from > 2.5<r peaks, constraining our analysis to 
22 haloes composed of > 3.2 x 10 5 particles and 29 haloes com- 
posed of > 1.2 x 10 5 particles for the n = —1 and n = —2.5 sim- 
ulation respectively. Each halo is associated with a velocity scale 



Aw = (GM h /R) 1/2 where R 3 = Afh/(47rp bg /- 3 /3), M h is 
the mass of the halo found using a linking length of £h and pbg 
is the background density. Two particles with phase-space coordi- 
nates (xi, vi) and (x2, V2) are linked if 



(xi - x 2 ) 2 (vi - Vg) 



< 1, 



(11) 



where £ s and b v are the dimensionless parameters defining the 
physical and velocity linking lengths respectively. We pass candi- 
date subhaloes through an unbinding routine that checks whether an 
object is self-bound and removes any unbound particles ( Springel 
|et al.|2 008 1. Only subhaloes with 20 particles or more are kept in 
the catalogue. This unbinding routine is a necessary but time con- 
suming process that removes unbound particles from subhalo can- 
didates and also eliminates spuriously linked particles. The frac- 
tion of unbound candidates in the initial 6DFOF catalogue varies 
slightly with b v and is « 0.05 and w 0.70 for the n = —1 and 
n = —2.5 simulation respectively. 

We set £ s = 0.10, thus limiting our search to re- 
gions within haloes with physical overdensities of p/pb g > 
1000, conversely we try several velocity linking lengths, b v — 
[0.0125,0.025,0.05,0.10]. For an isolated spherical overdensity, 
increasing b v amounts to placing a higher circular velocity cutoff in 
defining the boundary of a candidate subhalo. The mass associated 
with a phase-space peak would continuously increase as b v is in- 
creased were it not for our unbinding routine which effectively im- 
poses a tidal limit. In such a case, our technique is analogues to that 
used by Diemand et al. ( 2007 1. They found phase-space peaks using 
the 6DFOF algorithm and then estimated a tidal mass by assuming 
spherical symmetry and fitting an NFW profile plus background to 
circular velocity profiles of the peaks. Since our method does not 
impose a particular density profile, we can indirectly examine the 
validity of spherical symmetry asm —3 and any systematic bias 
this assumption introduces in the resulting subhalo mass function 
by varying b v . 

In Fig. [4] we show different representations of two large 
haloes, one from each of our simulations. The two haloes originated 
from ~ 3cr peaks and are composed of ~ 6 x 10 J particles. For 
each halo we show the phase-space density calculated using EnBiD 
( Sharma & Steinmetz 2006 1 and the subhalo candidates found us- 
ing two different velocity linking lengths. This figure demonstrates 
that there are significant substructure differences between the sim- 
ulations. Substructure at n = —2.5 appears more triaxial and less 
dense than at n — —1 and, overall, the n = — 1 halo looks very 
similar to a cluster or galaxy halo. At n = —1, substructure con- 
sists of distinct spherical overdensities with well defined bound- 
aries, whereas substructure is less evident at more negative indices, 
(see |Diemand et al.|2006| >. Of the 525 subhalo candidates found us- 
ing b v = 0.05, 464 are bound and contain a total of w 20% of the 
host halo's mass, though no single subhalo contains more than 5% 
of the halo's mass. In contrast, the n = —2.5 halo has 285 candi- 
dates, of which only 21 are bound. These bound subhaloes contain 
2% of the halo's mass. 

Comparing the impact of different velocity linking lengths b v , 
we find that at n = —2.5, using b v ^ 0.05 appears to spuriously 
link phase-space peaks into long, unbound, filamentary structures, 
which are subsequently eliminated from the catalogue by our un- 
binding routine. This is the case in the central region outlined by a 
solid black circle. However, not all peaks are linked into artificial 
filaments. Occasionally, peaks are close enough to one another in 
phase-space and are linked by more than a few particles. Two can- 
didate subhaloes, which are split at small b v , are linked at larger b v 
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Figure 3. The co-moving halo number density of mass M , dn(M, z)/d\r\ M denoted by filled points with lc error bars at n = — 1 (top) and n = —2.5 
(bottom). Curves are predictions from analytic fitting functions proposed by ST (solid) and PS (dashed). Going from right to left corresponds to going to 
progressively lower redshifts with right panel corresponding to the final redshift analyzed. 



into a bound triaxial object. An example of such is circled in a black 
dashed circle in Fig. [4] The dotted circle indicates a subhalo that is 
non-spherical regardless of b v used. These results visually demon- 
strate that the boundary of subhaloes appears less well defined as 

n — * — 3. 

5.2 Phase-space structure of haloes 

The upper panels of Fig. [4] provide the phase-space density as 
a function of position. Of course the phase-space distribution 
/(x, v), which provides a complete description of a dynamical 
system, is a function of six variables, making it cumbersome to 
deal with. As an alternative, we examine halo phase-space struc- 
ture using the volume distribution function «(/), which is the vol- 
ume of phase-space occupied by phase-space elements of density 
/ in an interval df ( |Arad et al.|2004| >. In Fig. [5] we plot v(f) (nor- 
malized so that J v(f)df = 1) along with the logarithmic slope, 
7„ = d\nv(f)/dln f for four haloes with > 3 x 10 5 particles. 
Haloes from the n = — 1 simulation have larger volumes with 
high phase-space density and span a greater range in / than the 
haloes from the n = —2.5 simulation. At n = — 1, 7„ decreases 
to « —2.5 with increasing / followed by a bump, and has a very 
similar form to that of ACDM galactic and cluster haloes examined 
by Sharma & Steinmetz (20061. The slope of haloes at n = —2.5 
is shallower at low / and plateaus around a slope of —2.5. Sev- 
eral authors have attempted to explain the form of 7„ for ACDM 



haloes using a toy model where haloes are modelled by a Hern- 
quist sphere that contains smaller Hernquist spheres representing 
subhaloes (e.g. lArad et al.]|2004| |Asc asibar & Binney 2003] and 
ISharma & Steinmetz|2 006 ). Sharma & Steinmetz ( 2006 ) found that 
the general shape of a ACDM galactic halo with substructure can 
be reproduced with these toy models and that the size of the bump 
appeared to be related to the amount of substructure. 

To examine whether the size of the bump in the slope is re- 
lated to the phase-space density contrast between substructure and 
the background, we smooth out substructure in the n = — 1 halo 
shown in Fig. [4] To smooth the halo we first calculate its morphol- 
ogy via the inertia tensor (see below) at different radii. Particles are 
then moved in a random direction through a small random angle 
on the surface of the ellipsoid determined by the mass distribution 
interior to their radii. This simple process leaves the radial density 
profile and overall morphology unchanged to within < 1% while 
reducing the density contrast of substructure, thereby decreasing 
the number of candidates and bound subhaloes by a factor of 2 — 3 
and 8 — 60 respectively. The result is shown in the right column 
of Fig. [5] We find that the size of the bump is influenced by both 
the amount and phase-space density contrast of substructure. The 
lack of a bump in the haloes from the n = —2.5 simulation shows 
that the contrast between substructure and the smooth background 
of a halo decreases as n — > —3. The differences in 7„ do not arise 
solely due to differences in halo morphology as the haloes in either 
simulation have similar morphologies. Thus, the logarithmic slope 
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Figure 4. Two haloes originating from 3cr peaks in the n = —1 (left) and n = — 2.5 simulation (right) composed of 6.2 X 10 5 and 5.5 X 10 5 particles 
respectively. The first row shows the phase-space density of the haloes using a logarithmic colour scale, where dark blue regions indicate high phase-space 
density.In the next two rows, particles are colour coded according to group number found with 6DFOF. Shown are subhalo candidates found with b v = 0.05 
and 0.0125 in the middle and bottom rows, respectively. Three regions are circled in the n = —2.5 simulation to highlight the variation in particles linked by 
the 6DFOF algorithm with different b v (see discussion in text). 
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Figure 5. The normalized volume distribution function v(f) (top) and logarithmic slope -~f v (bottom) for 4 haloes consisting of > 3 X 10 5 particles. The left 
and middle columns correspond to n = — 1 and n = —2.5 simulation respectively. The solid lines corresponds to the haloes shown in Fig. [4] the other line 
types are ordered in decreasing halo mass going from dashed, dotted to dashed-dotted. The right column correspond to a n = — 1 halo where the substructure 
has been smoothed with increasing smoothing going from dashed to dotted to dashed-dotted and the reference halo denoted by the solid curve. 



of v(f) can be a sensitive tool of the subhalo mass function and 
hence the differences seen in Fig. [5] are suggestive that there may 
be a spectral dependence in both the amplitude and the slope of 
subhalo mass function. 

5.3 Morphology 

To determine the morphology of a subhalo we follow [Dubinski 
|& Carlberg] l |1991| > and |Allgood et al.| ( [2006) and diagonalize the 
weighted moments of inertia tensor 

n 

The ellipsoidal distance between the subhalo's centre of mass and 
the nth particle is 

rl = x 2 n + (y„/q) 2 + {z n /sf, (13) 

where q and s are the intermediate-to-major and minor- to-major 
axis ratios respectively. The axis ratios are calculated after unbound 
particles are removed for subhaloes composed of > 100 particles. 

Figure [6] shows the distribution of subhaloes in terms of q and 
s. We see that subhaloes are generally more triaxial at n = —2.5 
compared to n = — 1 and the distribution of axis ratios is substan- 
tially broader. At both indices, the distribution remains relatively 
unchanged as b v is varied, save for the changes in the number of 



filamentary objects with q, s < 0.1. These filamentary objects are 
elongated subhaloes in the process of being tidally disrupted with 
loosely bound tidal tails. These filaments comprise < 1% of the 
subhalo population for b v ^ 0.05 at n = —1 and for b v ^ 0.025 
at n — —2.5, and increase to 20% for b v — 0.10 at both indices. 
This filamentary population decreases by roughly a factor of ~ 10 
and ~ 5 as b v is halved at n = —1 and n — —2.5 respectively, 
though the trend is not as clear cut at n — —2.5. The main reason 
for this dependence on b v is that tidal tails are more likely to be 
linked with the parent subhalo when using large b v . The decrease 
in the number of subhaloes as b v is increased at n = —2.5 indicates 
that, as n — > —3, local phase-space peaks become less well sepa- 
rated in phase-space and more likely to be embedded in a triaxial 
overdensity. The increasingly triaxial or filamentary nature makes 
assigning a mass to a phase-space peak non-trivial. Estimates of a 
tidal mass assuming spherical symmetry will become increasingly 
incorrect as the distribution of mass around a phase-space peak be- 
comes increasingly triaxial. These results are insensitive to the min- 
imum subhalo particle number cut applied. Again, the differences 
in morphology may be indicative of a break in the "universality" of 
the subhalo mass function. 
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Figure 6. Scatter plot of major and minor axis ratios, q and s, for subhaloes found in the n = —1 (top) and n = — 2.5 (bottom) simulations. Also shown 
are the projects of these distributions. The grey open circles, squares, triangles and upside down triangles and solid, dashed, dotted, and dashed-dotted lines 
correspond to b v = 0.0125, 0.025, 0.05 and 0.10 respectively. For clarity we only show a random subsample (10%) and plot contours for b v = 0.0125 and 
b v = 0.05. Filled points indicate the peaks the histograms and follow the same marker scheme as the scatter plot. The thin dashed line corresponds to q = s. 
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Figure 7. The cumulative subhalo mass function N(> Mf) (top) and logarithmic slope a (bottom) from all haloes used in our analysis. The line styles and 
marker scheme are the same as in Fig. [6] Also shown in the top panel by the solid thin black line is N(> M) oc M ~ The solid grey region outlined by 
black dashed lines in the bottom panel indicates the various slopes found in other studies with the solid horizontal line denoting a = —0.9. 



5.4 Mass function 

We define the dimensionless subhalo mass ratio Mf = 
Afsubhaio/Afhaio and show in Fig. [7] the cumulative subhalo 
mass function iV(> Mf) summed over all haloes in our 
analysis. We also show the logarithmic slope a(Mf) = 
dln(dN / din M f) / din Mf . Since dN/dlnMf is quite noisy, we 
calculate several logarithmic slopes at a given M f by varying the 
spacing used in the five point central difference estimate of the 
slope. Shown in Fig.[7]is the average slope along with the standard 
deviation. In both simulations, 7V(> Mf) is linear in the log-log 
plot provided we exclude low and high mass regions dominated 
by numerical effects. The flattening in the low mass region corre- 
sponds to subhaloes composed of fewer than 100 particles and is 
due to numerical softening effects. The high mass scale, where the 
distributions begins to steepen, scales roughly as 6 3 . This region re- 
sults from excluding the outer high velocity volumes of a candidate 
subhalo and the scaling can be understood as follows: v 2 oc M/ R, 
R oc A/ 1/3 , thus v 3 oc M. Hence, more massive the subhalo the 
larger b v must be to group particles out to the subhalo's tidal ra- 
dius. Since we pass subhalo candidates through an unbinding rou- 



tine, which does not take into account the background, the loosely 
unbound central regions of large subhaloes will be removed. 

In the n = — 1 simulation, the number of subhaloes is very 
weakly dependent on b v . The 6DFOF algorithm using b v — 0.10 
cannot separate subhaloes in the central region from the core and 
underestimates the number of subhaloes. The small changes in the 
number found with b v ^ 0.05 are due to the fact that the 6DFOF 
algorithm using these velocity linking lengths tends to only link the 
central regions of subhaloes with mass ratios of Mf > 0.01, un- 
derestimating their mass. In some cases, the algorithm only groups 
the loosely unbound centres of these subhaloes and they are subse- 
quently removed by our unbinding routine. However, this does not 
greatly affect the number since there are few subhaloes above the 
high mass limit imposed by b v . The fraction of mass bound in sub- 
haloes with 10~ 4 ^ Mf ^ 5 x 10 -3 is weakly dependent on b v , 
decreasing from 0.127 to 0.085 when going from b v — 0.0125 to 
b v = 0.10. This fraction varies from halo to halo by an average of 
« 10% for b v < 0.05 and increases to w 20% for b v = 0.10 with 
no strong dependence on peak height of the halo. The slope of the 
subhalo mass function is also unaffected by changes in b v and, ne- 
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glecting the regions dominated by numerical effects, is consistent 
with previous results, as indicated by the grey region in Fig. [7] Even 
between the catalogues found using b v = 0.10 and b v ^ 0.05, 
where half as many subhaloes are found and the fraction of fila- 
ments goes from 0.2 to < 0.01, the slope is only shallower for 
mass ratios of Mf < 2 x 10~ 4 . 

The subhalo distribution at n = —2.5 is dependent on b v . The 
number of subhaloes increases as b v decreases and appears to con- 
verge, though we note that few if any bound subhaloes are found us- 
ing b v < 0.0125. The increase in the number found is due to bound 
filamentary objects found at a given b v being broken into several 
smaller triaxial subhaloes at smaller b v , as shown in Fig. [4]&[6] 
For example, comparing the subhaloes found using b v — 0.025 to 
those found using b v — 0.0125, the number of filaments decreases 
from 1.7% to 0.3% of the population but the number of subhaloes 
increases by a factor of 1.6. The apparent convergence is due to 
there being a finite number of phase-space peaks with local physi- 
cal overdensities above > 1000. The fraction of mass in subhaloes 
between 10 -4 ^ Mf ^ 5 x 10 -3 shows a strong dependence on 
b v , decreasing from 0.035 to 0.0052 when going from b v = 0.0125 
to b v = 0.10. In general, the fraction of mass bound in subhaloes 
is a factor of ~ 3 — 20 smaller than at n — — 1. The relative halo 
to halo variation is also greater, with the fraction varying by 26% 
and 96% using b v = 0.0125 and b v — 0.10 respectively. The slope 
increases as b v increases and exhibits fluctuations as a function of 
Mf that are greater than those at n — — 1. Generally, a > —0.9 
within the region free of numerical effects. Given the changes in the 
number of subhaloes found, the changes in a could be dismissed 
as being purely numerical. However, similar changes in the sub- 
halo population are observed at n = — 1 going from b v = 0.10 to 
b v = 0.05 yet a is unaffected, whereas at n — —2.5 a increases by 
» 20% when going from b v = 0.05 to b v = 0.025 or b v = 0.025 
to b v — 0.0125. These changes in a, while being partly numeri- 
cal, are nonetheless indicative of the underlying physical issue of 
defining the boundary of a subhalo. 

To compare with previous work, we fit a power-law, 

N(> M f ) = AMf, (14) 

to the cumulative subhalo mass function, neglecting the low and 
high mass regions which are dominated by numerical effects. The 
mass scale above which subhaloes are artificially truncation due to 
b v is very clear at n — —1 but not as clear at n — —2.5, so we 
use the n = — 1 simulation for guidance in defining the range in 
Mf fitted. We plot the mean a as a function of the spectral index 
for different b v in Fig. [8] At n = — 1, a w —0.9 independent of 
the choice of velocity parameter in agreement with previous results 
from ACDM simulations. Examining the four largest haloes in the 
n = — 1 simulation, we find the halo to halo variation in the fitted 
index is 0.07. Unlike at n = —1, the index at n = —2.5 exhibits 
a strong dependence on b v , though it has similar a halo to halo 
variation of 0.06. We can reproduce a ~ —0.9 by choosing b v « 
0.0125, though in general a > —0.9. Increasing b v from 0.0125 
to 0.025, that is applying higher velocity cutoffs and then imposing 
a tidal limit, does flatten the slope by 22%. These results indicate 
that the slope of the subhalo mass function depends on n, becoming 
shallower but more sensitive to systematics asrn —3. 

To check for any redshift dependence, we examine the evolu- 
tion of a subset of our haloes over 0.8 and 0.5 e-foldings in expan- 
sion factors for 7i = —1 and n = —2.5 respectively. We also com- 
pare haloes between the two simulations at similar redshifts. How- 
ever, due to the redshift and mass dependence of v, we cannot at 
the same time continue to compare similar mass haloes originating 
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Figure 8. Spectral dependence of a. Filled stars, circles, squares 
and triangles indicate slopes from subhaloes found using b v = 
0.0125, 0.025, 0.05 and 0.10 respectively. Marker scheme is the same 
as in Fig. [6] The grey region and black lines indicate the same slopes as in 
Fig.^Jat n c ff = —1.8 to —2.2, which correspond to cluster scales down to 
galaxy scales in the CDM hierarchy. Error bars correspond to the average 
variation in the logarithmic slope seen in bottom panel of Fig. ^] within the 
region fitted by the power-law. 

from similar a peaks. We find that subhaloes do not have a strong 
systematic redshift dependence in either cosmology, thus the dif- 
ferences between the simulations remain unchanged. At n = — 1, 
a ~ —0.9 regardless of redshift and at n — —2.5 the value of 
a and its dependence on b v remains unchanged. It is interesting to 
note that a does vary randomly with redshift by w 0.05 and w 0.06 
for n = — 1 and n = —2.5 respectively. This variation is of similar 
size to the halo to halo variation seen, which for the subset studied 
spans a peak height range of 1 < v < 4.5 and 2.2 < v < 3.3 for 
n = — 1 and n = —2.5 respectively. 



6 DISCUSSION AND CONCLUSION 

One of the most intriguing results from the |Moore et al.| (1999| > 
study of CDM structure is similarity of substructure in galactic and 
cluster haloes. This result suggests that the subhalo mass distribu- 
tion is independent of parent halo mass, and possibly independent 
of shape of the primordial power spectrum. At first glance, this sim- 
ilarity may seem suprising since galaxy and cluster haloes differ by 
two orders of magnitude in mass. However, galaxies and clusters 
probe similar effective spectral indices, (n c ff ~ —2.2 as compared 
to riefr ~ —1-8), and it is the effective index that governs the struc- 
ture formation process via the relative formation times and merger 
rates of progenitors. 

To explore this issue, we have analyzed substructure in 
Einstein-de Sitter cosmologies at two widely separated indices, 
n — —1 and —2.5, which bracket the indices most often exam- 
ined in ACDM simulations. In the spirit of |Smith et al] {2003) 
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we have outlined a set of criteria to minimize the numerical ef- 
fects which plague simulations of n < —2 cosmologies. We chose 
n = —2.5 as it is the most negative index we can simulate that 
satisfies our criteria, albeit just barely. Another reason for choos- 
ing n = —2.5 is that it represents a transitional index between 
hierarchical formation and the singular, perhaps pathological, case 
of n — —3 where structures collapse simultaneously. We evolved 
our simulations until there were haloes large enough to analyze for 
substructure, while at the same time not being too significantly af- 
fected by finite volume effects. Substructure was analyzed with a 
number of tools: we examined the phase-space structure of haloes 
using EnBiD, searched for candidate subhaloes using a 6DFOF al- 
gorithm, and examined the morphology and mass distribution of 
bound subhaloes. All of these methods showed that substructure 
in an n = —2.5 cosmology is different from substructure in an 
n — —1 cosmology. For example, substructure leaves a feature 
in the volume distribution function whose strength depends on the 



number of subhaloes and the density contrast of substructure I Arad 



et al.||2004||Ascasibar & Binney 2005] and |Sharma & Stei nmetz 



2006). This feature was present in our n = — 1 simulation but was 
greatly reduced in our n — —2.5 simulation. 

The differences in the morphology of bound subhaloes be- 
tween the two cosmologies are striking and appear to be inde- 
pendent of redshift. At n = — 1, subhaloes were typically well- 
defined, spherical overdensities that would lend themselves to the 
type of analysis outlined by Diemand et al. ( 2007 ), where the tidal 
mass is determined by fitting spherically averaged density profiles 
to phase-space peaks. At n = —2.5, subhaloes were more triaxial 
and phase-space peaks were more irregular, making it difficult to 
define a boundary. The triaxial or filamentary nature of subhaloes, 
along with the large fraction of unbound candidates, indicates that 
substructure in an n = —2.5 cosmology is more susceptible to 
tidal disruption. These observations suggest that the methods that 
assume a spherical profile in order to estimate subhalo circular ve- 
locities and masses, such as those outlined in Diemand et al. ( 2006 
|2007[ l, will become increasingly inaccurate as n — > —3 and might 
give misleading results. 

The logarithmic slope of the subhalo mass distribution also 
showed a spectral dependence. At n = —1, we found a w —0.9 
independent of the the velocity linking length b v and is consistent 
with slopes from previous studies. At n — —2.5, we found a range 
of values —0.9 < a < —0.5 depending on our choice of b v with a 
tending towards less negative values with larger values of b v . The 
sensitivity of a to b v at n — —2.5 is due to the underlying un- 
certainty in determining the boundary of a subhalo. Our preferred 
value is a — —0.76 ± 0.07, found using b v = 0.025, for a num- 
ber of reasons: subhaloes were found over several decades in Mf 
below the numerical mass limit imposed by b v and above the mass 
scales dominated by numerical softening effects; highly filamen- 
tary subhaloes made up « 1% of the population; and the number of 
subhaloes appeared to be converging. We also find that the logarith- 
mic slope of the mass function shows no dependence on redshift, 
though a does vary by ~ 5 — 10% in time for a given halo. This 
random variation in a with redshift is of similar size halo to halo 
variation, which shows little dependence on the peak height of the 
halo. 

Our preferred slope at n = —2.5 appears to agree with that 
of the subsubhalo mass function in ~ 10 9 Mq subhaloes studied 
by Springel et al. (2008), where n c ff « —2.4. However, we note 
that the subsubhalo mass functions in Springel et al. (2008]) are 
quite noisy with a great deal of scatter in the observed logarithmic 
slope, so this agreement should be treated with caution. This value 



is not entirely consistent with that found in Gao et al. ( 2003} and 
|Diemand et al.| {2006), who both found results consistent with pre- 
vious studies of much larger scales. Using SUBFIND, |Gao et ak] 
{2005} examined substructure in a single high redshift subgalac- 
tic halo at n t fs < —2.5 originating from an extremely rare m 5er 
peak, whereas we examined the logarithmic slope averaged over a 
number of haloes, albeit they examined their halo over a larger red- 
shift range than we did. At their earliest epoch, this rare peak corre- 
sponded to a 2 x 10 5 Mq halo composed of ~ 2 x 10 5 particles with 
n c ff w —2.7. Though their halo sampled a more negative index then 
our simulation, it was an extremely rare peak that only had w 35 
subhaloes over a sing le decade in Mf from 10~ 4 < Mf < 10" 3 . 
We also note that the logarithmic slope of this rare peak appears to 
vary by ~ 0.1 with redshift, though there is no systematic depen- 
dence on redshift, similar to the variation we observed. Di emand] 
et al. (2006) also examined a single rare halo, though at an effec- 
tive index slightly closer to —3 and at a much higher resolution than 
|Gao et ah] {2005 ) study. They noted that the subhalo mass function 
was sensitive to the parameters used to find candidates, unlike the 
mass function of larger haloes where n e tf > —2, in agreement with 
our observations, although they used SKID and we used a 6DFOF. 
Our results, when combined with those of previous studies 
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mic slope, though relativity constant for a given host, is dependent 
on the host's mass and its associated spectral index. This slope is 
constant with a value of ~ —0.9 so long as n > —2 but tends to 
larger values and possibly with increasing scatter for n < —2. This 
transition in behaviour might be due to the structure formation pro- 
cess. Provided that this process is "sufficiently" hierarchical, that 
is haloes virialize before being accreted or merging, the subhalo 
mass function is independent of n. As n — + —3, haloes do not 
fully virialize before being accreted and consequently the subhalo 
distribution is influenced by the properties of the power spectrum. 
This dependence manifests itself in the boundary of a subhalo be- 
coming increasingly ill defined as n — > —3; phase-space peaks 
become more indistinct and irregular, resulting in greater variation 
in the subhalo mass function while also tending to flatten it. Thus, 
extrapolating the subhalo mass function at galactic scales to pre- 
dict the number of subhaloes at the bottom of the CDM hierarchy 
is questionable. It is likely that at n = — 3 subhaloes become gen- 
erally impossible to distinguish and, in essence, substructure be- 
comes negligible. 

The distribution and internal properties of subhaloes have im- 
portant ramifications for dark matter detectors. Because the 7-ray 
flux varies as the local dark matter density squared, indirect detec- 
tors, such as GLAST, will see a stronger signal if there is a great 
deal of dense substructure. Numerous groups have calculated the 
boost in flux due to substructure (e.g. |Strigari et al.||2007| |Pieri| 
et al. 2008, and Kuhlen et al. 2008) but there are two issues with 
these estimates. The primary issue is the scale dependence of the 
subhalo mass function. Estimates assume the subhalo mass func- 
tion observed in 10 12 Mq haloes at subhalo masses of > 10 5 Mq 
applies to the subsubhalo population and extends to the bottom of 
the CDM hierarchy. However, it is likely that the subsubhalo mass 
function of ~ 10 s Mq subhaloes will have a form similar to that 
observed at n = —2.5 for b v = 0.025 and that the distribution will 
continue to flatten as one proceeds down the CDM hierarchy. The 
secondary issue is the use of field halo mass-concentration relations 
to convert subhalo masses to densities. At small scales, their use is 
questionable since subhaloes become increasingly triaxial and less 
distinct as n — > —3. However, this could be accounted for in the 
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concentration-mass relation such as the one proposed by Eke et al. 
( 200 1 1 which attempts to account for the possible effects of n c ff . Ul- 
timately, as a consequence of these issues, previous results should 
generally be considered optimistic upper limits. In light of these re- 
sults we do not attempt to predict the 7-ray background here and 
will address this critical issue in a paper in preparation. 
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